#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkPoint.h"
#include "itkImageRegionIteratorWithIndex.h"
#include <fstream>
#include <vector>
#include <list>
#include <string>

// Constants
#define DIMENSION 3
#define PHI 1.017221967897852 // constant to generate icosahedron
#define THETA 0.553574358897045
#define R1 1.902113032590307
#define R2 1.095445115010332
#define R3 1.129775730952839
#define PI 3.141592653589793
// Macro : a, b are both 3D vectors, vecprod( a, b ) compute their dot product
#define vecprod( a, b ) \
	( a.GetX() * b.GetX() + a.GetY() * b.GetY() + a.GetZ() * b.GetZ() )
#define get_x( i ) ( tetra->Vert[i].GetX() )
#define get_y( i ) ( tetra->Vert[i].GetY() )
#define get_z( i ) ( tetra->Vert[i].GetZ() )
#define set_x( i, x ) ( tetra->Vert[i].SetX( ( x ) ) )
#define set_y( i, y ) ( tetra->Vert[i].SetY( ( y ) ) )
#define set_z( i, z ) ( tetra->Vert[i].SetZ( ( z ) ) )
#define vec( i ) ( tetra->Vert[i] )
#define sq( a ) ( ( a ) * ( a ) )
// typedef
typedef unsigned char										PixelType ;
typedef itk::Image< PixelType, DIMENSION >					ImageType ;
typedef itk::Point< double, DIMENSION >						PointType ;

class NodeElem 
{
public:
	NodeElem()
	{
		this->m_x = 0.0 ;
		this->m_y = 0.0 ;
		this->m_z = 0.0 ;
		this->m_density = 0.0 ;
		this->m_index = 0 ;
	}
	NodeElem( float X, float Y, float Z, float Density, int Index )
	{
		this->m_x = X ;
		this->m_y = Y ;
		this->m_z = Z ;
		this->m_density = Density ;
		this->m_index = Index ;
	}
	NodeElem( char* X, char* Y, char* Z, char* Density, int Index )
	{
		this->m_x = atof( X ) ;
		this->m_y = atof( Y ) ;
		this->m_z = atof( Z ) ;
		this->m_density = atof( Density ) ;
		this->m_index = Index ;
	}
	double GetX()
	{
		return m_x ;
	}
	double GetY()
	{
		return m_y ;
	}
	double GetZ()
	{
		return m_z ;
	}
	double GetDensity()
	{
		return m_density ;
	}
	void SetX( double X )
	{
		m_x = X ;
	}
	void SetY( double Y )
	{
		m_y = Y ;
	}
	void SetZ( double Z )
	{
		m_z = Z ;
	}
	void ConvertMesh2RAW( int sizeY, double spacingY )
	{
		double MeshY = this->m_y ;
		this->m_y = this->m_x ;
		this->m_x = sizeY * spacingY - MeshY ;
	}
private:
	double m_x, m_y, m_z, m_density ;
	int m_index ;
};

struct SrcTetraElem
{
	NodeElem Vert[4] ;
	double Cen_X, Cen_Y, Cen_Z, a, b, c, Dist[4] ;
	int index_a, index_b, index_c ;
};

struct DetectIcoElem
{
	int sumvox, vesvox ;
	double Cen_X, Cen_Y, Cen_Z ;
};

void ComputeEllipsoidParam( SrcTetraElem* tetra )
{
	// initialize
	tetra->Cen_X = tetra->Cen_Y = tetra->Cen_Z
		= tetra->a = tetra->b = tetra->c = 0.0 ;
	tetra->index_a = tetra->index_b = tetra->index_c = 0 ;

	// Compute the centroid
	for ( int i = 0 ; i < 4 ; i++ )
	{
		tetra->Cen_X += get_x(i) ;
		tetra->Cen_Y += get_y(i) ;
		tetra->Cen_Z += get_z(i) ;
	}
	tetra->Cen_X /= 4.0 ;
	tetra->Cen_Y /= 4.0 ;
	tetra->Cen_Z /= 4.0 ;

	// Compute a ( the longest semi-principal axis )
	for ( int i = 0 ; i < 4 ; i++ )
	{
		set_x( i, get_x( i ) - tetra->Cen_X ) ;
		set_y( i, get_y( i ) - tetra->Cen_Y ) ;
		set_z( i, get_z( i ) - tetra->Cen_Z ) ;
		tetra->Dist[i] = sqrt( vecprod( vec(i), vec(i) ) ) ;
		tetra->a = tetra->a > tetra->Dist[i] ? tetra->a : tetra->Dist[i] ;
		tetra->index_a = tetra->a > tetra->Dist[i] ? tetra->index_a : i ;
	}

	NodeElem inda = vec( tetra->index_a ) ;

	// Compute the norm of the other three vector right to a
	for ( int i = 0; i < 4 ; i++ )
	{
		if ( i != tetra->index_a )
		{
			NodeElem indi = vec(i) ;
			set_x( i, get_x( i ) - vecprod( indi, inda ) /
				vecprod( inda, inda ) * inda.GetX() ) ;
			set_y( i, get_y( i ) - vecprod( indi, inda ) /
				vecprod( inda, inda ) * inda.GetY() ) ;
			set_z( i, get_z( i ) - vecprod( indi, inda ) /
				vecprod( inda, inda ) * inda.GetZ() ) ;
			tetra->Dist[i] = sqrt( vecprod( vec( i ), vec( i ) ) ) ;
			tetra->b = tetra->b > tetra->Dist[i] ? tetra->b : tetra->Dist[i] ;
			tetra->index_b = tetra->index_b > tetra->Dist[i] ? tetra->index_b : i ;
		}
	}

	NodeElem indb = vec( tetra->index_b ) ;

	// Compute the norm of the other two vector to a and b
	for ( int i = 0; i < 4 ; i++ )
	{
		if ( i != tetra->index_a && i != tetra->index_b )
		{
			NodeElem indi = vec(i) ;
			set_x( i, get_x( i ) - vecprod( indi, indb ) /
				vecprod( indb, indb ) * indb.GetX() ) ;
			set_y( i, get_y( i ) - vecprod( indi, indb ) /
				vecprod( indb, indb ) * indb.GetY() ) ;
			set_z( i, get_z( i ) - vecprod( indi, indb ) /
				vecprod( indb, indb ) * indb.GetZ() ) ;
			tetra->Dist[i] = sqrt( vecprod( vec( i ), vec( i ) ) ) ;
			tetra->c = tetra->Dist[i] ;
			tetra->index_c = i ;
			break ;
		}
	}
}

void IcosahedronInit( std::list<DetectIcoElem*>& icosahedronquad )
{
	// top
	// 1, ( 0, 0, 1 )
	DetectIcoElem* vertex = new DetectIcoElem ;
	vertex-> Cen_X = 0.0 ;
	vertex-> Cen_Y = 0.0 ;
	vertex-> Cen_Z = 1.0 ;
	icosahedronquad.push_back( vertex ) ;
	// mid of top and positive pentagon
	// 2, ( cos(0.0*pi), sin(0.0*pi), cos(theta)*r1/r1 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = 1.0 / R1 ;
	vertex-> Cen_Y = 0.0 ;
	vertex-> Cen_Z = cos( THETA ) * R1 / R1 ;
	icosahedronquad.push_back( vertex ) ;
	// 3, ( cos(0.4*pi), sin(0.4*pi), cos(theta)*r1/r1 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.4 * PI ) / R1 ;
	vertex-> Cen_Y = sin( 0.4 * PI ) / R1 ;
	vertex-> Cen_Z = cos( THETA ) ;
	icosahedronquad.push_back( vertex ) ;
	// 4, ( cos(0.8*pi), sin(0.8*pi), cos(theta)*r1/r1 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.8 * PI ) / R1 ;
	vertex-> Cen_Y = sin( 0.8 * PI ) / R1 ;
	vertex-> Cen_Z = cos( THETA ) ;
	icosahedronquad.push_back( vertex ) ;
	// 5, ( cos(1.2*pi), sin(1.2*pi), cos(theta)*r1/r1 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.2 * PI ) / R1 ;
	vertex-> Cen_Y = sin( 1.2 * PI ) / R1 ;
	vertex-> Cen_Z =  cos( THETA ) ;
	icosahedronquad.push_back( vertex ) ;
	// 6, ( cos(1.6*pi), sin(1.6*pi), cos(theta)*r1/r1 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.6 * PI ) / R1 ;
	vertex-> Cen_Y = sin( 1.6 * PI ) / R1 ;
	vertex-> Cen_Z = cos( THETA ) ;
	icosahedronquad.push_back( vertex ) ;
	// mid of top pentagon
	// 7, (cos(0.2*pi), sin(0.2*pi), cos(phi)/r3 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.2 * PI ) / R3 ;
	vertex-> Cen_Y = sin( 0.2 * PI ) / R3 ;
	vertex-> Cen_Z = cos( PHI ) / R3 ;
	icosahedronquad.push_back( vertex ) ;
	// 8, ( cos(0.6*pi), sin(0.6*pi), cos(phi)/r3 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.6 * PI ) / R3 ;
	vertex-> Cen_Y = sin( 0.6 * PI ) / R3 ;
	vertex-> Cen_Z = cos( PHI ) / R3 ;
	icosahedronquad.push_back( vertex ) ;
	// 9, ( cos(1.0*pi), sin(1.0*pi), cos(phi)/r3 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = -1.0 / R3 ;
	vertex-> Cen_Y =  0.0 ;
	vertex-> Cen_Z = cos( PHI ) / R3 ;
	icosahedronquad.push_back( vertex ) ;
	// 10,( cos(1.4*pi), sin(1.4*pi), cos(phi)/r3  )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.4 * PI ) / R3 ;
	vertex-> Cen_Y = sin( 1.4 * PI ) / R3 ;
	vertex-> Cen_Z = cos( PHI ) / R3 ;
	icosahedronquad.push_back( vertex ) ;
	// 11,( cos(1.8*pi), sin(1.8*pi), cos(phi)/r3 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.8 * PI ) / R3 ;
	vertex-> Cen_Y = sin( 1.8 * PI ) / R3 ;
	vertex-> Cen_Z = cos( PHI ) / R3 ;
	icosahedronquad.push_back( vertex ) ;
	// positive pentagon
	// 12,( cos(0.0*pi), sin(0.0*pi), cos(2*theta)/r2 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = 1.0 / R2 ;
	vertex-> Cen_Y = 0.0 ;
	vertex-> Cen_Z = cos( 2 * THETA ) / R2 ;
	icosahedronquad.push_back( vertex ) ;
	// 13,( cos(0.4*pi), sin(0.4*pi), cos(2*theta)/r2 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.4 * PI / R2 ) ;
	vertex-> Cen_Y = sin( 0.4 * PI / R2 ) ;
	vertex-> Cen_Z = cos( 2 * THETA ) / R2 ;
	icosahedronquad.push_back( vertex ) ;
	// 14,(cos(0.8*pi), sin(0.8*pi), cos(2*theta)/r2 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.8 * PI ) / R2 ;
	vertex-> Cen_Y = sin( 0.8 * PI ) / R2 ;
	vertex-> Cen_Z = cos( 2 * THETA ) / R2 ;
	icosahedronquad.push_back( vertex ) ;
	// 15,( cos(1.2*pi), sin(1.2*pi), cos(2*theta)/r2 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.2 * PI ) / R2 ;
	vertex-> Cen_Y = sin( 1.2 * PI ) / R2 ;
	vertex-> Cen_Z = cos( 2 * THETA ) / R2 ;
	icosahedronquad.push_back( vertex ) ;
	// 16,( cos(1.6*pi), sin(1.6*pi), cos(2*theta)/r2 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.6 * PI ) / R2 ;
	vertex-> Cen_Y = sin( 1.6 * PI ) / R2 ;
	vertex-> Cen_Z = cos( 2 * THETA ) / R2 ;
	icosahedronquad.push_back( vertex ) ;
	// equitor
	// 17,( cos(0.1*pi), sin(0.1*pi), 0 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.1 * PI ) ;
	vertex-> Cen_Y = sin( 0.1 * PI ) ;
	vertex-> Cen_Z = 0.0 ;
	icosahedronquad.push_back( vertex ) ;
	// 18,( cos(0.3*pi), sin(0.3*pi), 0 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.3 * PI ) ;
	vertex-> Cen_Y = sin( 0.3 * PI ) ;
	vertex-> Cen_Z = 0.0 ;
	icosahedronquad.push_back( vertex ) ;
	// 19,( cos(0.5*pi), sin(0.5*pi), 0 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.5 * PI ) ;
	vertex-> Cen_Y = sin( 0.5 * PI ) ;
	vertex-> Cen_Z = 0.0 ;
	icosahedronquad.push_back( vertex ) ;
	// 20,( cos(0.7*pi), sin(0.7*pi), 0 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.7 * PI ) ;
	vertex-> Cen_Y = sin( 0.7 * PI ) ;
	vertex-> Cen_Z = 0.0 ;
	icosahedronquad.push_back( vertex ) ;
	// 21,( cos(0.9*pi), sin(0.9*pi), 0 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.9 * PI ) ;
	vertex-> Cen_Y = sin( 0.9 * PI ) ;
	vertex-> Cen_Z = 0.0 ;
	icosahedronquad.push_back( vertex ) ;
	// 22,( cos(1.1*pi), sin(1.1*pi), 0 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.1 * PI ) ;
	vertex-> Cen_Y = sin( 1.1 * PI ) ;
	vertex-> Cen_Z = 0.0 ;
	icosahedronquad.push_back( vertex ) ;
	// 23,( cos(1.3*pi), sin(1.3*pi), 0 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.3 * PI ) ;
	vertex-> Cen_Y = sin( 1.3 * PI ) ;
	vertex-> Cen_Z = 0.0 ;
	icosahedronquad.push_back( vertex ) ;
	// 24,( cos(1.5*pi), sin(1.5*pi), 0 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.1 * PI ) ;
	vertex-> Cen_Y = sin( 0.1 * PI ) ;
	vertex-> Cen_Z = 0.0 ;
	icosahedronquad.push_back( vertex ) ;
	// 25,( cos(1.7*pi), sin(1.7*pi), 0 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.7 * PI ) ;
	vertex-> Cen_Y = sin( 1.7 * PI ) ;
	vertex-> Cen_Z = 0.0 ;
	icosahedronquad.push_back( vertex ) ;
	// 26,( cos(1.9*pi), sin(1.9*pi), 0 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.9 * PI ) ;
	vertex-> Cen_Y = sin( 1.9 * PI ) ;
	vertex-> Cen_Z = 0.0 ;
	icosahedronquad.push_back( vertex ) ;
	// negative pentagon
	// 27,( cos(0.2*pi), sin(0.2*pi), -cos(2*theta)/r2 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.2 * PI ) / R2 ;
	vertex-> Cen_Y = sin( 0.2 * PI ) / R2 ;
	vertex-> Cen_Z = -cos( 2 * THETA ) / R2 ;
	icosahedronquad.push_back( vertex ) ;
	// 28,( cos(0.6*pi), sin(0.6*pi), -cos(2*theta)/r2 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.6 * PI ) / R2 ;
	vertex-> Cen_Y = sin( 0.6 * PI ) / R2 ;
	vertex-> Cen_Z = -cos( 2 * THETA ) / R2 ;
	icosahedronquad.push_back( vertex ) ;
	// 29,( cos(1.0*pi), sin(1.0*pi), -cos(2*theta)/r2 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.0 * PI ) / R2 ;
	vertex-> Cen_Y = sin( 1.0 * PI ) / R2 ;
	vertex-> Cen_Z = -cos( 2 * THETA ) / R2 ;
	icosahedronquad.push_back( vertex ) ;
	// 30,( cos(1.4*pi), sin(1.4*pi), -cos(2*theta)/r2 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.4 * PI ) / R2 ;
	vertex-> Cen_Y = sin( 1.4 * PI ) / R2 ;
	vertex-> Cen_Z = -cos( 2 * THETA ) / R2 ;
	icosahedronquad.push_back( vertex ) ;
	// 31,( cos(1.8*pi), sin(1.8*pi), -cos(2*theta)/r2 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.8 * PI ) / R2 ;
	vertex-> Cen_Y = sin( 1.8 * PI ) / R2 ;
	vertex-> Cen_Z = -cos( 2 * THETA ) / R2 ;
	icosahedronquad.push_back( vertex ) ;
	// mid of bottom pentagon
	// 32,( cos(0.0*pi), sin(0.0*pi), -cos(phi)/r3 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = 1.0 / R3 ;
	vertex-> Cen_Y = 0.0 ;
	vertex-> Cen_Z = -cos( PHI ) / R3 ;
	icosahedronquad.push_back( vertex ) ;
	// 33,( cos(0.4*pi), sin(0.4*pi), -cos(phi)/r3 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.4 * PI ) / R3 ;
	vertex-> Cen_Y = sin( 0.4 * PI ) / R3 ;
	vertex-> Cen_Z = -cos( PHI ) / R3 ;
	icosahedronquad.push_back( vertex ) ;
	// 34,( cos(0.8*pi), sin(0.8*pi), -cos(phi)/r3 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.8 * PI ) / R3 ;
	vertex-> Cen_Y = sin( 0.8 * PI ) / R3 ;
	vertex-> Cen_Z = -cos( PHI ) / R3 ;
	icosahedronquad.push_back( vertex ) ;
	// 35,( cos(1.2*pi), sin(1.2*pi), -cos(phi)/r3 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.2 * PI ) / R3 ;
	vertex-> Cen_Y = sin( 1.2 * PI ) / R3 ;
	vertex-> Cen_Z = -cos( PHI ) / R3 ;
	icosahedronquad.push_back( vertex ) ;
	// 36,( cos(1.6*pi), sin(1.6*pi), -cos(phi)/r3 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.6 * PI ) / R3 ;
	vertex-> Cen_Y = sin( 1.6 * PI ) / R3 ;
	vertex-> Cen_Z = -cos( PHI ) / R3 ;
	icosahedronquad.push_back( vertex ) ;
	// mid of bottom and negative pentagon
	// 37,( cos(0.2*pi), sin(0.2*pi), -cos(theta)*r1/r1 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.2 * PI ) / R1 ;
	vertex-> Cen_Y = sin( 0.2 * PI ) / R1 ;
	vertex-> Cen_Z = -cos( THETA ) ;
	icosahedronquad.push_back( vertex ) ;
	// 38,( cos(0.6*pi), sin(0.6*pi), -cos(theta)*r1/r1 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 0.6 * PI ) / R1 ;
	vertex-> Cen_Y = sin( 0.6 * PI ) / R1 ;
	vertex-> Cen_Z = -cos( THETA ) ;
	icosahedronquad.push_back( vertex ) ;
	// 39,( cos(1.0*pi), sin(1.0*pi), -cos(theta)*r1/r1 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = -1.0 / R1 ;
	vertex-> Cen_Y =  0.0 ;
	vertex-> Cen_Z = -cos( THETA ) ;
	icosahedronquad.push_back( vertex ) ;
	// 40,( cos(1.4*pi), sin(1.4*pi), -cos(theta)*r1/r1 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.4 * PI ) / R1 ;
	vertex-> Cen_Y = sin( 1.4 * PI ) / R1 ;
	vertex-> Cen_Z = -cos( THETA ) ;
	icosahedronquad.push_back( vertex ) ;
	// 41,( cos(1.8*pi), sin(1.8*pi), -cos(theta)*r1/r1 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = cos( 1.8 * PI ) / R1 ;
	vertex-> Cen_Y = sin( 1.8 * PI ) / R1 ;
	vertex-> Cen_Z = -cos( THETA ) ;
	icosahedronquad.push_back( vertex ) ;
	// bottom
	// 42, ( 0, 0, -1 )
	vertex = new DetectIcoElem ;
	vertex-> Cen_X = 0.0 ;
	vertex-> Cen_Y = 0.0 ;
	vertex-> Cen_Z = -1.0 ;
	icosahedronquad.push_back( vertex ) ;
}

void IcosahedronVertCoord( std::list<DetectIcoElem*>& icosahedronquad,
	double radius )
{
	std::list<DetectIcoElem*>::iterator icosahedronIter ;
	for ( icosahedronIter = icosahedronquad.begin() ;
		icosahedronIter != icosahedronquad.end() ; icosahedronIter++ )
	{
		( *icosahedronIter)->Cen_X *= radius ;
		( *icosahedronIter)->Cen_Y *= radius ;
		( *icosahedronIter)->Cen_Z *= radius ;
	}
}

void CountVesselVoxels( std::list<DetectIcoElem*>::iterator Iter,
	ImageType::Pointer volume, double Cen_X, double Cen_Y,	double Cen_Z,
	double radius, int lowerlim, int upperlim )
{
	( *Iter )->sumvox = 0 ;
	( *Iter )->vesvox = 0 ;
	PointType detectcenter, regionStart, regionEnd ;
	detectcenter[0] = Cen_X + ( *Iter )->Cen_X ;
	detectcenter[1] = Cen_Y + ( *Iter )->Cen_Y ;
	detectcenter[2] = Cen_Z + ( *Iter )->Cen_Z ;
	regionStart = detectcenter - radius ;
	regionEnd = detectcenter + radius ;
	typedef itk::NearestNeighborInterpolateImageFunction<
		ImageType, double > InterpolatorType ;
	InterpolatorType::Pointer intpltor = InterpolatorType::New() ;
	intpltor->SetInputImage( volume ) ;
	
	typedef itk::ImageRegionIteratorWithIndex< ImageType > IteratorType ;
	ImageType::RegionType region ;
	ImageType::RegionType::IndexType idxStart, idxEnd ;
	intpltor->ConvertPointToNearestIndex( regionStart, idxStart ) ;
	intpltor->ConvertPointToNearestIndex( regionEnd, idxEnd ) ;
	ImageType::RegionType::SizeType volsize =
		volume->GetBufferedRegion().GetSize() ;
	if ( idxEnd[0] < 0 || idxEnd[1] < 0 || idxEnd[2] < 0 ||
		idxStart[0] > volsize[0] || idxStart[1] > volsize[1] ||
		idxStart[2] > volsize[2] )
	{
		return ;
	}
	idxStart[0] = idxStart[0] < 0 ? 0 : idxStart[0] ;
	idxStart[1] = idxStart[1] < 0 ? 0 : idxStart[1] ;
	idxStart[2] = idxStart[2] < 0 ? 0 : idxStart[2] ;
	idxEnd[0] = idxEnd[0] > volsize[0] ? volsize[0] : idxEnd[0] ;
	idxEnd[1] = idxEnd[1] > volsize[1] ? volsize[1] : idxEnd[1] ;
	idxEnd[2] = idxEnd[2] > volsize[2] ? volsize[2] : idxEnd[2] ;

	ImageType::RegionType::SizeType regionSize ;
	regionSize[0] = idxEnd[0] - idxStart[0] + 1 ;
	regionSize[1] = idxEnd[1] - idxStart[1] + 1 ;
	regionSize[2] = idxEnd[2] - idxStart[2] + 1 ;
	region.SetSize( regionSize ) ;
	region.SetIndex( idxStart ) ;
	IteratorType regionIter( volume, region ) ;
	ImageType::SpacingType spacing = volume->GetSpacing() ;
	for ( regionIter.GoToBegin(); !regionIter.IsAtEnd(); ++regionIter )
	{
		PointType pnt ;
		ImageType::RegionType::IndexType idx ;
		pnt[0] = regionIter.GetIndex()[0] * spacing[0] ;
		pnt[1] = regionIter.GetIndex()[1] * spacing[1] ;
		pnt[2] = regionIter.GetIndex()[2] * spacing[2] ;
		double dst = sq( pnt[0] - detectcenter[0] ) +
			sq( pnt[1] - detectcenter[1] ) + sq( pnt[2] - detectcenter[2] ) ;
		if ( dst < sq( radius ) )
		{
<<<<<<< .mine
			intpltor->ConvertPointToNearestIndex( pnt, idx ) ;
=======
			intpltor->ConvertPointToNearestIndex( pnt, idx ) ;
			volume->SetPixel( idx, 240 ) ;
>>>>>>> .r11
			( *Iter )->sumvox++ ;
			if (regionIter.Get() >= lowerlim &&
				regionIter.Get() <= upperlim )
			{
				( *Iter )->vesvox++ ;
<<<<<<< .mine
				volume->SetPixel( idx, 50 ) ; // overlay
=======
				volume->SetPixel( idx, 230 ) ;
>>>>>>> .r11
			}
			else
        	{
	    		volume->SetPixel( idx, 30 ) ;// small ball
        	}
		}
	}
}

void ComputeVesselIntensity( std::list<DetectIcoElem*>& icosahedronquad,
	ImageType::Pointer volume, std::list<SrcTetraElem*> SrcTetraList,
	double radius, int lowerlim, int upperlim )
{
	double Cen_X, Cen_Y, Cen_Z, vertDen ;
	Cen_X = Cen_Y = Cen_Z = vertDen = 0.0 ;
	std::list<SrcTetraElem*>::iterator SrcTetraListIter = SrcTetraList.begin() ;
	for ( SrcTetraListIter = SrcTetraList.begin() ;
		SrcTetraListIter != SrcTetraList.end() ; SrcTetraListIter ++  )
	{
		for ( int i = 0 ; i < 4 ; i++ )
		{
			if ( ( *SrcTetraListIter )->Vert[i].GetDensity() > vertDen )
			{
				vertDen = ( *SrcTetraListIter )->Vert[i].GetDensity() ;
				Cen_X = ( *SrcTetraListIter )->Vert[i].GetX() +
					( *SrcTetraListIter )->Cen_X ;
				Cen_Y = ( *SrcTetraListIter )->Vert[i].GetY() +
					( *SrcTetraListIter )->Cen_Y ;
				Cen_Z = ( *SrcTetraListIter )->Vert[i].GetZ() +
					( *SrcTetraListIter )->Cen_Z ;
			}
		}
	}
	std::list<DetectIcoElem*>::iterator icosahedronquadIter
		= icosahedronquad.begin() ;
	for ( icosahedronquadIter = icosahedronquad.begin() ;
		icosahedronquadIter != icosahedronquad.end() ; icosahedronquadIter++ )
	{
		CountVesselVoxels( icosahedronquadIter, volume, Cen_X, Cen_Y, Cen_Z,
			radius,	lowerlim, upperlim ) ;
	}
}

void main( int argc, char* argv[] )
{
	typedef itk::ImageFileReader< ImageType >					ReaderType ;
	typedef itk::ImageFileWriter< ImageType >					WriterType ;

	ImageType::Pointer rawVolume ;
	ReaderType::Pointer reader = ReaderType::New() ;
	reader->SetFileName( argv[1] ) ;
	rawVolume = reader->GetOutput() ;
	std::cout << "Reading volume data from file " << argv[1] << " ..." << std::endl ;
	try
	{
		reader->Update() ;
	}
	catch( itk::ExceptionObject & excep )
	{
		std::cerr << "Exception caught !" << std::endl ;
		std::cerr << excep << std::endl ;
	}
	std::cout << "volume data have been read successfully!" << std::endl ;

	ImageType::RegionType rawRegion = rawVolume->GetLargestPossibleRegion() ;
	ImageType::SpacingType rawSpacing = rawVolume->GetSpacing() ;
	ImageType::PointType rawOrigin = rawVolume->GetOrigin() ;
	ImageType::RegionType::IndexType rawStart = rawRegion.GetIndex() ;
	ImageType::SizeType rawSize = rawRegion.GetSize() ;
	double minSp = rawSpacing[0] ; // minimum spacing
	for ( int i = 1; i < DIMENSION; i++ )
	{
		minSp = minSp < rawSpacing[i] ? minSp : rawSpacing[i] ;
	}

	std::ifstream datfile( argv[2] ) ;
	if ( datfile.bad() )
	{
		return ;
	}
	std::cout << "Reading mesh data from file " << argv[2] << " ..." << std::endl ;
	// read the tetrahedra from file

	char title[50] ;
	datfile.getline( title, 50, '\n' ) ;
	char category[50] ;
	datfile.getline( category, 50, '\n' ) ;
	char discrib[50] ;
	datfile >> discrib ;
	datfile.getline( discrib, 50, '\n' ) ;
	int N ;
	sscanf( discrib, " N=%d", &N ) ;
	char data[50] ;
	float X, Y, Z, Density ;
	X = Y = Z = Density = 0.0f ;
	
	std::vector< NodeElem > NodeVector ;
	std::list< NodeElem > SrcNodeList ;
	NodeVector.reserve( N ) ;
	for ( int i = 1 ; i <= N ; i++ )
	{
		datfile.getline( data, 50, '\n' ) ;
		sscanf( data, "%f %f %f %f", &X, &Y, &Z, &Density ) ;
		NodeElem elem = NodeElem( X, Y, Z, Density, i ) ;
		elem.ConvertMesh2RAW( rawSize[1], rawSpacing[1] ) ;
		NodeVector.push_back( elem ) ;
		if ( Density != 0.0f )
		{
			SrcNodeList.push_back( elem ) ;
		}
	}

	std::cout << "Computing ellipsoids parameters and rendering..." << std::endl ;
	char tetraNodeIndex[5] ;
	std::list<SrcTetraElem*> SrcTetraList ;
	datfile >> tetraNodeIndex ;
	while ( tetraNodeIndex[0] != 'Z' )
	{
		// Compute the parameters of the ellipsoids
		SrcTetraElem* tetra = new SrcTetraElem ;
		tetra->Vert[0] = NodeVector[ atoi( tetraNodeIndex ) - 1 ] ;
		datfile >> tetraNodeIndex ;
		tetra->Vert[1] = NodeVector[ atoi( tetraNodeIndex ) - 1 ] ;
		datfile >> tetraNodeIndex ;
		tetra->Vert[2] = NodeVector[ atoi( tetraNodeIndex ) - 1 ] ;
		datfile >> tetraNodeIndex ;
		tetra->Vert[3] = NodeVector[ atoi( tetraNodeIndex ) - 1 ] ;
		SrcTetraList.push_back( tetra ) ;
		datfile >> tetraNodeIndex ;
		ComputeEllipsoidParam( tetra ) ;
	}

	// editing center
	std::list<SrcTetraElem*>::iterator SrcTetraListIter = SrcTetraList.begin() ;
	( *SrcTetraListIter )->Cen_X +=  6.72;
	( *SrcTetraListIter )->Cen_Y += -17.28;
	( *SrcTetraListIter )->Cen_Z += -1.32;
	SrcTetraListIter ++ ;
	( *SrcTetraListIter )->Cen_X +=  6.72;
	( *SrcTetraListIter )->Cen_Y += -17.28;
	( *SrcTetraListIter )->Cen_Z += -1.32;
	// ...

	// rendering the tetrahedra in the volume data
	for ( SrcTetraListIter = SrcTetraList.begin() ; SrcTetraListIter != SrcTetraList.end() ; SrcTetraListIter ++  )
	{
		typedef itk::NearestNeighborInterpolateImageFunction<
			ImageType, double > InterpolatorType ;
		InterpolatorType::Pointer intpltor = InterpolatorType::New() ;
		intpltor->SetInputImage( rawVolume ) ;
		PointType point ;
		point.Fill( 0.0 ) ;
		ImageType::RegionType::IndexType idx ;
		idx.Fill( 0 ) ;
		
		int gra, grb, grc ;
		SrcTetraElem* tetra = *SrcTetraListIter ;
		for ( gra = 0 ;	gra < (int)( 2 * tetra->a / minSp ) ; gra++ )
		{
			point[0] = tetra->Cen_X - tetra->a + ( ( double )gra ) *
				minSp * get_x( tetra->index_a ) / tetra->a ;
			point[1] = tetra->Cen_Y - tetra->b + ( ( double )gra ) *
				minSp * get_y( tetra->index_a ) / tetra->a ;
			point[2] = tetra->Cen_Z - tetra->c + ( ( double )gra ) *
				minSp * get_z( tetra->index_a ) / tetra->a ;
			for ( grb = 0 ; grb < (int)( 2 * tetra->b / minSp ) ; grb++ )
			{
				for ( grc = 0 ; grc < (int)( 2 * tetra->c / minSp ) ; grc++ )
				{
					if ( sq( ( ( ( ( double )gra ) * minSp - tetra->a ) / tetra->a ) )
						+ sq( ( ( ( ( double )grb ) * minSp - tetra->b ) / tetra->b ) )
						+ sq( ( ( ( ( double )grc ) * minSp - tetra->c ) / tetra->c ) )
						<= 1.0 )
					{
						intpltor->ConvertPointToNearestIndex( point, idx ) ;
<<<<<<< .mine
						rawVolume->SetPixel( idx, 10 ) ; //stem cell
=======
						rawVolume->SetPixel( idx, 255 ) ;
>>>>>>> .r11
					}
					point[0] += minSp *	get_x( tetra->index_c ) / tetra->c ;
					point[1] += minSp *	get_y( tetra->index_c ) / tetra->c ;
					point[2] += minSp *	get_z( tetra->index_c ) / tetra->c ;
				}
				point[0] -= grc * minSp * get_x( tetra->index_c ) / tetra->c ;
				point[1] -= grc * minSp * get_y( tetra->index_c ) / tetra->c ;
				point[2] -= grc * minSp * get_z( tetra->index_c ) / tetra->c ;
				point[0] += minSp * get_x( tetra->index_b ) / tetra->b ;
				point[1] += minSp * get_y( tetra->index_b ) / tetra->b ;
				point[2] += minSp * get_z( tetra->index_b ) / tetra->b ;
			}
		}
	}
	std::cout << "volume data rendering finished!" << std::endl ;

	// Output rendered volume
	WriterType::Pointer writer = WriterType::New() ;
	writer->SetInput( rawVolume ) ;
	std::string filename( argv[1] ) ;
	filename.insert( filename.length() - 4, "_rendered" ) ;
	writer->SetFileName( filename ) ;
	try
	{
		writer->Update() ;
	}
	catch( itk::ExceptionObject & excep )
	{
		std::cerr << "Exception caught !" << std::endl ;
		std::cerr << excep << std::endl ;
	}
	std::cout << "Rendered Volume written to file : " << filename << std::endl ;

	// Compute vessel density
	int lowerlim = atoi( argv[3] ) ; // the lower limit of vessel gray intensity in the volume data
	int upperlim = atoi( argv[4] ) ; // the upper limit of vessel gray intensity in the volume data
	double r_Diffuse = atof( argv[5] ) ; // the radius of the distribution of detectors
	double r_Detect = atof( argv[6] ) ; // the radius of the detecting area
	std::list<DetectIcoElem*> icosahedron ;
	std::list<DetectIcoElem*>::iterator icoIter ;
	IcosahedronInit( icosahedron ) ;
	IcosahedronVertCoord( icosahedron, r_Diffuse ) ;
	ComputeVesselIntensity( icosahedron, rawVolume, SrcTetraList, r_Detect,
		lowerlim, upperlim ) ;
	double vesselDensity = 0.0, sumVol = 0.0, vesselVol = 0.0 ;
	for ( icoIter = icosahedron.begin() ; icoIter != icosahedron.end();
		++icoIter )
	{
		sumVol += ( *icoIter )->sumvox ;
		vesselVol += ( *icoIter )->vesvox ;
	}
	vesselDensity = vesselVol / sumVol * 100 ;

	try
	{
		writer->Update() ;
	}
	catch( itk::ExceptionObject & excep )
	{
		std::cerr << "Exception caught !" << std::endl ;
		std::cerr << excep << std::endl ;
	}


	try
	{
		writer->Update() ;
	}
	catch( itk::ExceptionObject & excep )
	{
		std::cerr << "Exception caught !" << std::endl ;
		std::cerr << excep << std::endl ;
	}

	std::cout << "Vessel density around the stem cells is "
		<< vesselDensity << "%.( " << vesselVol << " voxels in " << sumVol << " voxels)" << std::endl ;
}